!-----------
       module param
       implicit none

       integer :: lx,ly,nb,nh,nn,mm
       real*8 :: beta,aprib,bprob,aprob,dprob
       real*8, external :: drand1,ran
       real*8 :: enrg1, enrg2,amag2,ususc
       integer, dimension(:), allocatable :: spin,frstspinop,lastspinop,vertexlist,opstring
       integer, dimension(:,:), allocatable :: bsites

       end

!-----------
       program sse_heisenberg
       use param
       implicit none

       integer :: i,j,nbins,msteps,isteps,iseed

       namelist /lattice/ lx,ly
       namelist /choice/ beta,nbins,msteps,isteps,iseed

!       open(unit=10,file='input.d',form='formatted',status='unknown')

!      lattice
       lx=0
       ly=0
       read(5,lattice)

       beta=2.d0
       nbins=100
       msteps=10000
       isteps=10000
       iseed=235716731
       read(5,choice)

       nh=0
       nn=lx*ly
       mm=max(4,nn/4)
       allocate(spin(nn),frstspinop(nn),lastspinop(nn),opstring(0:mm-1),vertexlist(0:4*mm-1))
       opstring=0

       enrg1=0.d0
       enrg2=0.d0
       amag2=0.d0 
       ususc=0.d0 

       call rand_init(iseed)

       call makelattice()

       call initconfig()

       aprob=0.5d0*beta*nb
       dprob=1.d0/(0.5d0*beta*nb)

       do i=1,isteps
          call diagonalupdate()
          call linkvertices()
          call loopupdate()

          call adjustcutoff(i)
       end do

       do j=1,nbins
          do i=1,msteps
             call diagonalupdate()
             call linkvertices()
             call loopupdate()
             call measure()
          end do
          call writeresults(msteps)
       end do

       stop
       end

!-----------
       subroutine makelattice()
       use param
       implicit none

       integer :: s,x1,x2,y1,y2

       if(ly==1)then
         nb=lx
         allocate(bsites(2,nb))
         bsites=0
         do x1=0,lx-1
            s=x1+1
            x2=mod(x1+1,lx)
            bsites(1,s)=s
            bsites(2,s)=x2+1
!            write(*,*) x1,bsites(1,s),bsites(2,s)
         end do

       else if(ly==2)then

         nb=nn+lx
         allocate(bsites(2,nb))
         bsites=0
         do y1=0,ly-1
            do x1=0,lx-1
               s=1+x1+y1*lx
               x2=mod(x1+1,lx)
               y2=y1
               bsites(1,s)=s
               bsites(2,s)=1+x2+y2*lx
               if(y1==0)then
                 x2=x1
                 y2=mod(y1+1,ly)
                 bsites(1,s+nn)=s
                 bsites(2,s+nn)=1+x2+y2*lx       
               end if
!               write(*,*) x1,y1,bsites(1,s),bsites(2,s),bsites(1,s+nn),bsites(2,s+nn)
            end do
         end do

       else

         nb=2*nn
         allocate(bsites(2,nb))
         bsites=0
         do y1=0,ly-1
            do x1=0,lx-1
               s=1+x1+y1*lx
               x2=mod(x1+1,lx)
               y2=y1
               bsites(1,s)=s
               bsites(2,s)=1+x2+y2*lx
               x2=x1
               y2=mod(y1+1,ly)
               bsites(1,s+nn)=s
               bsites(2,s+nn)=1+x2+y2*lx 
!               write(*,*) x1,y1,bsites(1,s),bsites(2,s),bsites(1,s+nn),bsites(2,s+nn)
            end do
         end do
       end if

       continue
       end

!-----------------------!
       subroutine initconfig()
       use param
       implicit none

       integer :: i

       do i=1,nn
          spin(i)=(-1)**(mod(i-1,lx)+(i-1)/lx)
       end do
!       write(*,*) 'spin', spin

       continue
       end

!---------------------------!
       subroutine diagonalupdate()
       use param
       implicit none

       integer :: i,b,op

!       write(*,*) 'before diag'
!       write(*,*) 'spin',spin
!       write(*,*) 'opstring',opstring

       do i=0,mm-1
          op=opstring(i)
          if(op==0)then       
            b=int(drand1()*nb)+1
!            b=int(ran()*nb)+1
!            write(*,*) 'inside',i,b
            if(spin(bsites(1,b))/=spin(bsites(2,b)))then
              if(drand1()*(mm-nh)<=aprob)then
!              if(ran()*(mm-nh)<=aprob)then
                opstring(i)=2*b
                nh=nh+1 
              end if
            end if
          else if(mod(op,2)==0)then        
            if(drand1()<=dprob*(mm-nh+1))then
!            if(ran()<=dprob*(mm-nh+1))then
              opstring(i)=0
              nh=nh-1
            end if
          else
            b=op/2
            spin(bsites(1,b))=-spin(bsites(1,b))
            spin(bsites(2,b))=-spin(bsites(2,b))
          end if
       end do
!       write(*,*) 'after diag'
!       write(*,*) 'spin',spin
!       write(*,*) 'opstring',opstring

       continue
       end

!-------------------------!
       subroutine linkvertices()
       use param
       implicit none

       integer :: b,op,s1,s2,v0,v1,v2

       frstspinop=-1
       lastspinop=-1

!       write(*,*) 'before link'
!       write(*,*) 'vertex',vertexlist
!       write(*,*) 'frst',frstspinop
!       write(*,*) 'last',lastspinop

       do v0=0,4*mm-1,4
          op=opstring(v0/4)
          if(op/=0)then
            b=op/2
            s1=bsites(1,b)
            s2=bsites(2,b)
            v1=lastspinop(s1)
            v2=lastspinop(s2)
!            write(*,*) 'test',v0,op,b,s1,s2,v1,v2

            if(v1/=-1)then
              vertexlist(v1)=v0
              vertexlist(v0)=v1
            else
              frstspinop(s1)=v0
            end if

            if(v2/=-1)then
              vertexlist(v2)=v0+1
              vertexlist(v0+1)=v2
            else
              frstspinop(s2)=v0+1
            end if
            lastspinop(s1)=v0+2
            lastspinop(s2)=v0+3
          else
            vertexlist(v0:v0+3)=-1
          end if

!          write(*,*) 'inside link',v0
!          write(*,*) 'vertex',vertexlist
!          write(*,*) 'frst',frstspinop
!          write(*,*) 'last',lastspinop
       end do

       do s1=1,nn
          v1=frstspinop(s1)
          if(v1/=-1)then
            v2=lastspinop(s1)
            vertexlist(v2)=v1
            vertexlist(v1)=v2
          end if
       end do

!       write(*,*) 'after link'
!       write(*,*) 'vertex',vertexlist
!       write(*,*) 'frst',frstspinop
!       write(*,*) 'last',lastspinop

       continue
       end

!-----------------------!
       subroutine loopupdate()
       use param
       implicit none

       integer :: i,v0,v1,v2

!       write(*,*) 'before loop'
!       write(*,*) 'vertex',vertexlist
!       write(*,*) 'spin',spin

       do v0=0,4*mm-1,2       
          if(vertexlist(v0)<0) cycle
          v1=v0
          if(drand1()<0.5d0)then
            do 
               opstring(v1/4)=ieor(opstring(v1/4),1)
               vertexlist(v1)=-2
               v2=ieor(v1,1)
               v1=vertexlist(v2)
               vertexlist(v2)=-2
               if(v1==v0) exit
            end do
          else
            do 
               vertexlist(v1)=-1
               v2=ieor(v1,1)
               v1=vertexlist(v2)
               vertexlist(v2)=-1
               if(v1==v0) exit
            end do
          end if
       end do

       do i=1,nn
          if(frstspinop(i)/=-1)then
            if(vertexlist(frstspinop(i))==-2) spin(i)=-spin(i)
          else
            if(drand1()<0.5) spin(i)=-spin(i)
          end if
       end do

!       write(*,*) 'after loop'
!       write(*,*) 'vertex',vertexlist
!       write(*,*) 'spin',spin

       continue
       end

!-----------------------------!
 subroutine adjustcutoff(step)
!--------------------------------------------------------------------!
! Increases the cut-off mm of the expansion in case the current mm is
! less than nh+nh/3 (nh being the number of H operators in the string).
! The arrays of size mm must be deallocated and allocated again.
!--------------------------------------------------------------------!

 use configuration; implicit none

 integer, allocatable :: stringcopy(:)
 integer :: mmnew,step

 mmnew=nh+nh/3
 if (mmnew<=mm) return

 allocate(stringcopy(0:mm-1))
 stringcopy(:)=opstring(:)

 deallocate(opstring)
 allocate(opstring(0:mmnew-1))
 opstring(0:mm-1)=stringcopy(:)
 opstring(mm:mmnew-1)=0
 deallocate(stringcopy)

 mm=mmnew
 deallocate (vertexlist)
 allocate(vertexlist(0:4*mm-1))

       open(unit=10,file='cut.dat',position='append')
       write(10,*)step,mm
       close(10)

       continue
       end

!--------------------!
       subroutine measure()
       use param
       implicit none

       integer :: i,b,op,s1,s2,am
       real(8) :: am2

       am=0
       do i=1,nn
          am=am+spin(i)*(-1)**(mod(i-1,lx)+(i-1)/lx)
       end do
       am=am/2
       am2=0.d0
       do i=0,mm-1
          op=opstring(i)
          if(mod(op,2)==1)then
            b=op/2
            s1=bsites(1,b)
            s2=bsites(2,b)
            spin(s1)=-spin(s1)
            spin(s2)=-spin(s2)
            am=am+2*spin(s1)*(-1)**(mod(s1-1,lx)+(s1-1)/lx)
          end if
          am2=am2+dfloat(am)**2
       end do
       am2=am2/dble(mm)

       enrg1=enrg1+dble(nh)
       enrg2=enrg2+dble(nh)**2
       amag2=amag2+am2
       ususc=ususc+dble(sum(spin)/2)**2

       continue
       end
!----------------------!

!----------------------!
       real(8) function ran()
       implicit none

       real(8)    :: dmu64
       integer(8) :: ran64,mul64,add64
       common/bran64/dmu64,ran64,mul64,add64

       ran64=ran64*mul64+add64
       ran=0.5d0+dmu64*dble(ran64)

       end function ran

!---------------------!
       subroutine initran(w)
       implicit none

       integer(8) :: irmax
       integer(4) :: w,nb,b

       real(8)    :: dmu64
       integer(8) :: ran64,mul64,add64
       common/bran64/dmu64,ran64,mul64,add64

       irmax=2_8**31
       irmax=2*(irmax**2-1)+1
       mul64=2862933555777941757_8
       add64=1013904243
       dmu64=0.5d0/dble(irmax)
       open(10,file='seed.in',status= 'old')
       read(10,*)ran64
       close(10)

       if (w.ne.0) then
          open(10,file='seed.in',status='unknown')
          write(10,*)abs((ran64*mul64)/5+5265361)
          close(10)
       endif

       end subroutine initran
